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Molecular Dynamic (MD) approach is applied to study the converging cylindrical shock 
waves in a dense Lennard- Jones (LJ) fluid. MD method is based on tracking of the atom motions 
and hence it has an fundamental advantages over hydrodynamic methods which assumes shocks 
as a structureless discontinuity and requires an equation of state. Due to the small thickness of 
shock fronts in liquid the two million particles is enough to simulate propagation of a cylindrical 
shocks in close detail. 

We investigate stability of converging shocks with different perturbation modes and its mix- 
ture. It was shown that in a case of relatively large initial ripples the Mach stems are formed. Su- 
personic jets generated by interaction of reflected shocks in downstream flow are observed. We 
also study the Richtmyer-Meshkov (RM) instability of an interface between two Lennard- Jones 
liquids of different mass densities. Surprisingly, mode 3 ripples grow very slow in comparison 
with higher mode numbers and growth rate of a higher mode decay slower. 

1. Introduction 

Molecular Dynamic simulations involving the converging shock stability and the Richtmyer- 
Meshkov fluid instability at atomic scale are presented for the first time. Understanding of 
these instabilities is very important in inertial fusion science. So far, the capacity of the fluid 
numerical methods to compute fluid flow at small lengths calls in question, especially in case 
of converging shocks OX 

Last 40 years MD method was successfully applied to many problems in different fields 
such as statistical and chemical physics. With increasing of the computer power, simulations 
of many millions of atoms in 3D space are now feasible. We believe that nowadays is time to 
apply MD method to fluid hydrodynamic problems. 

MD approach is based on tracking of the atom motions and hence it can provide us whole 
information of a system behavior including phase transition in the matter, diffusion, thermo- 
conductivity and other nonequilibrium processes far from thermodynamic steady state. Due 
to direct computation of the atom movements and the atom-atom interactions MD method has 
an fundamental advantages over numerical hydrodynamic methods which assumes shocks as a 
structureless discontinuity and requires an equation of state. Moreover, the different hydrody- 
namic codes show significant deviations in the simulation flow and a higher mesh resolution can 
not improve the situation [2], but the various MD programs give the same simulation results. 

Fluid simulations of the Richtmyer-Meshkov instability by using numerical hydrodynamic 
methods have so far failed to provide quantitative agreement with experiments due to grid in- 
duced numerical instability. For this reason, it is necessary to apply an alternate method, namely 
MD approach as a very promising method of hydrodynamic simulation. It is the purpose of this 
work to carry out such numerical experiments. Simulations were performed by own full vector- 
ized MD code - Molecular Dynamics Solver (11 seconds per one time step for 2 21 LJ atoms, 
and 92 ns per a pair of atoms on the NEC SX-5 supercomputer ). 



Fig. 1 MD simulation of perturbed shock front at m = 8 and v p = 1.25 . Brightness cor- 
responds to atom density, each pixel is occupied by 20 atoms roughly, (left, t = 38.6 ): First 
generation of Mach stems appears at very early stage of converging SW. (center, t = 50.5 ): 
Phase inversion of the perturbations occurs. (right, t = 96.5 ): Reflected shock has almost un- 
perturbed symmetrical front and 8 supersonic jets run away from the target surface. Animation 
is located on www.ile.osaka-u.ac.jp/research/TSI/Vasilii/ , file nd-pb08.gif, size 5.1 MB. 



2. Simulation technique and results 

In our MD simulation the atoms of a liquid interact via modified Lennard- Jones short-range 
pair potential with cut-off and smoothing at r cut = 2.5a J3]|: 



LJ (r) = 4e (a/r) 12 - (a/rf - a 2 (r 2 - r 2 ) 2 - a 3 (r 2 - r 2 ) 3 , if < r < r, 



cut J 



where r = 2 1 / 6 cr , and a , e are the usual LJ parameter, and a 2 = —3.5289 x 10~ 3 , a 3 = 
5.75868 x 1CT 4 . We use these parameters and the atomic mass m a /48 as reduced molecular 
dynamic units. Here and after all quantities are in MD units (mdu). For argon atoms a = 



3.405 A, e/k B = 119.8 K, time mdu = a^m a /48e = 3.113 x 10~ 13 s, velocity mdu = 1094 
m/s, and the unit of pressure is 0.0419 GPa . 

At the beginning atoms are placed into rectangular MD cell inside cylindrical domain sur- 
rounded by a piston as a cylinder wall. The piston is simulated by an external potential ~ 
[r — R((f), t)} 2 and its position depends on angles in case of a perturbed boundary as well as 
on time to generate shock waves. The total number of atoms in simulation is 2097148 . The 
computational MD cell has dimensions L x x L y x L z , where L x = L y = 549.7 and the 
thickness of MD cell is L z = 21.8 , and periodical boundary conditions are imposed on the 
system along the cylinder axis (z-axis). 

At the preparation stage the cylinder radius is fixed R = 197.25 and all atoms have the 
same atomic weight m a = 48 . The Langevin thermostat is applied to prepare initial cold LJ 
liquid at the thermodynamic equilibrium with given temperature T = 0.72 and number density 
n = 0.79 . In case of a perturbed cylinder boundary the initial position of the piston is defined 
as R{4>) = Rq[1 + 5sm(m(j))] where m is the mode number, and 5 = 0.05 is the relative 
perturbation, A = 25 Ro is the initial perturbation amplitude. 

In the end of the preparation stage an equilibrium state with uniform mass density is reached. 
At the beginning of the simulation stage in order to simulate two different materials the atoms 
outside of the interface change instantly its mass to heavy one nib = 16m a and its velocities to 
Vb = v a /A to hold uniform temperature and pressure. Interaction between light-heavy atoms 




and heavy-heavy atoms is described by the same LJ potential (1) without any changing of the 
parameters. At t — the piston starts to move with the constant velocity v p and at t — 16 the 
piston is removed from the system. We estimate the sound velocity in light material as c = 0.7 
and measured shock Mach number is M = 3.3 for piston velocity v p — 1 . 




Fig. 3 MD simulation of RM instability for m = 8 and v p — 1 . Brightness corresponds 
to atom density, only heavy material is shown, (left, t = 68.3 ): Shock front hits material 
interface. Reflected wave is rarefaction wave, (center, t = 139.5 ): Spikes and bubbles begin 
to grow after hitting of the interface by SW reflected from origin (see Fig. 5, second triangle) 
(right, t = 267.2 ): Mixing zone mounts almost the asymptotic stage. Animation is located on 
www.ile.osaka-u.ac.jp/research/TSWasilii/ , nd-rm08.gif, size 8.2 MB. 

Typical snapshots of MD system in the case of the perturbed converging shock are shown in 
Figure 1 . It is found that for enough large initial perturbation the curved shock front generates 
Mach stems at very early stage of converging. Then the first generation of Mach stems entails 
the secondary Mach stems and so on. In fact polygons converge to the center instead of smooth 
curved shock front. There are m -side polygons for mode m or 2m -sides at phase inversion 
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Fig. 4 Flow velocity map at heavy-light interface in case of the RM instability for m = 8 . 
Dark arrows correspond to heavy material, light arrows - light material. Corresponding density 
maps are shown in Fig. 3. (left, t = 68.3 ): SW-interface interaction. Each arrow is an averaging 
of 320 atom velocities roughly. Frame size is 70 mdu. (right, t = 139.5 ): Flow velocity map 
around spike. Each arrow represents the mean velocity about 100 atoms. Frame size is 40 mdu. 



stage. Also the Mach stems generate intricate structures of reflected shocks in downstream flow 
which result in supersonic jets running away from the target boundary. We suppose these jets 
can amplify instability of ablation front in the real ICF experiments. Figure 2 shows amplitudes 
of the measured shock front perturbation as a function of a mean shock radius and the perturba- 
tion growth rate for different mode numbers. The simulation results show that the perturbations 
grow if the mode number m < 5 @|. 

Figures 3 and 4 represent simulation results of the RM instability where SW propagates from 
heavy to light liquids. The typical flow velocity fields are presented in Fig. 4. When the shock 
front hits the interface, the transmitted shock and rarefaction wave are produced and go away 
from the material interface (Figs. 3,A:left). It is clearly seen that the vorticities are allocated 
around the material interface and they cause a spike to be formed (Fig.4:right). ¥ig.3:right 
shows the RM instability at late simulation time when the total width of the mixing zone goes 
into the asymptotic stage of the growth. 

The time history of the RM instability is shown in Figure 5. We measure the thickness of 
the mixing zone A as a distance between the maximum and minimum of interface positions 
along the radius. For comparison with the symmetrical case we draw mean position (circles) of 
the cylindrical interface m = as a function of time. It should be mentioned that small atomic 
perturbations with m ^> 1 even in a pure symmetrical case are unavoidable due to atomistic 
fluctuations and they begin grow as well as artificial ripples (see A in Fig. 5 left). 

There are three key time points in the RM instability evolution. The first point (1st triangle 
in the left graph of Fig. 5) corresponds to a contact of the shock front with a bottom of the ripple. 
After that the maximum and minimum are closing in time. The phase inversion occurs at the 
minimum of A , between 1st and 2nd triangles. For the higher mode number the inversion 
takes place at earlier time than that for the lower mode number. The shock reflected at the 
origin reaches the interface at the second key time point and it accelerates the growth of the 
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Fig. 5 Time history of the RM instability. Rb and R s are maximum and minimum of 
interface location. At the late time, after 2nd triangle, they correspond to the radial position of 
the bubbles and spikes, (left): Circles denote interface position in symmetrical case without any 
ripples. Trajectory of the cylindrical shock front Rsw,o is shown because of perturbed shock 
has a close location on average. Triangles correspond to key points of SW-interface interaction 
(see text), (right): Growth rates of perturbations vs. time for different mode number. 



interface perturbations. In the case of mode 5 the coincidence of the phase inversion with the 
second point occurs accidentally, the phase inversion of the interface for m = 3 takes place 
after the secondary reflected shock pushes the interface outward at 3rd key time point, t ~ 180 . 

Figure 5:right indicates nonlinear evolution of the growth rates for the different modes. Due 
to numerical differentiation of the thickness of the mixing zone A = A(t) , the growth rate 
dA/dt is noisy and we are forced to estimate asymptotical behavior of the growth rate from 
the differentiation of A(t) asymptotic form. In the case of mode 8 the perturbations grow as 
A 8 ~ t 45 , and A 5 ~ t 3 , A 3 ~ ln(t) . Hence the growth rates decay dA 8 /dt ~ r°- 55 , and 
dA 5 /dt ~ t~ - 7 , and dA 3 /dt ~ t^ 1 correspondingly. 

3. Conclusion 

Molecular dynamics simulation of cylindrical shocks have been performed for the first time. 
We have demonstrated that MD method is able to catch hydrodynamic instabilities in close 
details on atomic scale. 

It is shown that for a converging cylindrical shock M 3.3 , the compression of Argon 
liquid is 2.2 at the collapse point, the temperature is around 12000 K, pressure is « 50 GPa, 
converging time 20 ps. The growth rate of a ripples strongly depends on the perturbation 
modes and slowly decreases in time. Perturbations grow faster for higher mode in the case of 
the RM instability. Mode 3 ripples grow very slow in comparison with higher mode numbers 
due to the phase inversion of the interface occurs at very late time. 
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